A simplified particulate model for coarse-grained hemodynamics simulations 
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Human blood flow is a multi-scale problem: in first approximation, blood is a dense suspension 
of plasma and deformable red cells. Physiological vessel diameters range from about one to thou- 
sands of cell radii. Current computational models either involve a homogeneous fluid and cannot 
track particulate effects or describe a relatively small number of cells with high resolution, but are 
incapable to reach relevant time and length scales. Our approach is to simplify much further than 
existing particulate models. We combine well established methods from other areas of physics in 
order to find the essential ingredients for a minimalist description that still recovers hemorheology. 
These ingredients are a lattice Boltzmann method describing rigid particle suspensions to account 
for hydrodynamic long range interactions and — in order to describe the more complex short-range 
behavior of cells — anisotropic model potentials known from molecular dynamics simulations. Paying 
detailedness, we achieve an efficient and scalable implementation which is crucial for our ultimate 
goal: establishing a link between the collective behavior of millions of cells and the macroscopic 
properties of blood in realistic flow situations. In this paper we present our model and demonstrate 
its applicability to conditions typical for the microvasculature. 

PACS numbers: 47.11. Qr, 87.19.U-, 83.10.Rs, 87.19. rh 



I. INTRODUCTION 

Human blood is not a homogeneous substance but can 
be approximated as a suspension of deformable red blood 
cells (RBCs, also called erythrocytes) in a Newtonian liq- 
uid, the blood plasma. Under physiological conditions 
the volume concentration of RBCs typically amounts to 
40 % to 50 % in larger blood vessels. We neglect the 
other constituents like leukocytes and thrombocytes in 
this work due to their far lower volume concentrations [l| . 
In the absence of external stresses, erythrocytes assume 
the shape of biconcave discs with a diameter of approxi- 
mately 8/j,m [3|. Their main biological task is the trans- 
port of oxygen in the body, but due to the high volume 
fraction they also strongly affect the rheology of blood 
and its clotting behavior [3| . An understanding of these 
effects is necessary in order to study and to cure patho- 
logically deviating phenomena in the body and to design 
microfluidic devices for improved blood analysis. In both 
cases, blood often has to be studied within complex ge- 
ometries that elude an analytical description. However, 
also the computational treatment of blood is demanding. 
On large scales like in arteries with diameters in the order 
of millimeters, blood can be modeled as a continuous and 
even Newtonian fluid [j|. Even then, the computational 
effort and the complexity of the model can be significant 
if realistic geometries show features which stretch over 
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different length scales. 

For modeling flow in the microvascular network, there 
is need for a description that accounts for the presence 
of discrete cells [3 • Recently, models of deformable cells 
were presented amongst others by Dupin et al. [6|, in 
the group of Gompper |7|, and by Wu and Aidun @. 
Here, the cell membrane is simplified to a deformable 
mesh and coupled to a mesoscopic simulation method 
for the plasma like multi-particle collision dynamics |7J or 
lattice Boltzmann (LB) [0, l3|. However, mainly because 
of the high resolution that is necessary for the elaborate 
description of the cells, these models are computationally 
too demanding for the application to considerably larger 
3D systems. 

Our motivation is to bridge the scales that are ac- 
cessible with both classes of existing models by an in- 
termediate approach: we keep the particulate nature of 
blood, but try to find a minimal description of each sin- 
gle cell. We thus deliberately simplify much further than 
the authors of the particulate models cited above in order 
to gain the potential for a computationally efficient and 
scalable implementation. Resorting to well-established 
methods from other areas of physics we explore the in- 
gredients necessary to recover the rheological behavior of 
blood. It is not our motivation to account for sub-cell 
effects in more than a coarse-grained way. In this work 
we aim at the description and validation of such a model 
while presenting possible applications in the range from 
approximately 10 yum upward. Our ultimate goal is to de- 
velop a quantitative method that allows to study the flow 
in realistic geometries but also to link bulk properties, for 
example the apparent viscosity, to phenomena at the level 
of single erythrocytes. In the case of cell deformation and 



aggregation in plane shear flow, this hnk has been estab- 
hshed ah'eady in experiment and theory (9|. Numerical 
simulations enable us to extend this knowledge to the 
case of arbitrary geometries and time-dependent flows. 
Further microscopic properties of interest are the align- 
ment of cells or local changes of the cell concentration. 
The improved understanding of the dynamic behavior of 
blood might be used for the optimization of macroscopic 
simulation methods. Only a computationally efficient de- 
scription allows the accumulation of firm statistical data 
that is necessary for this task. 

The main idea of our model is to distinguish between 
the long-range hydrodynamic coupling of cells and the 
short-range interactions that are related to the complex 
mechanics, electrostatics, and the chemistry of the mem- 
branes. Long-range hydrodynamic interactions are con- 
sidered by means of the LB method [13 • This mesoscopic 
simulation method allows a relatively easy implementa- 
tion of complex boundary conditions which are needed for 
the simulation of realistic geometries. Moreover, an ef- 
ficient parallelization is straightforward which even with 
our simplified model is crucial for the accumulation of 
statistically relevant data or for the description of real- 
istic systems like branching vessels. Research on a par- 
allel and efficient implementation of the LB method for 
the simulation of flow in s par se vessel networks was pub- 
lished by various authors [ll| - [l3| . Consequently, the LB 
method was applied to blood flow already in earlier works 
though they differ from our approach in the accuracy 
with which cells are resolved. For example, Boyd et al. 
[J| model blood either as a Newtonian fluid or a homo- 
geneous fluid with a shear rate-dependent viscosity. Fur- 
ther, the studies by Dupin et al. [6|, Wu and Aidun [8|, 
Migliorini et al. [1J|, and by Sun and Munn [15| involve 
LB solvers, but describe each RBC as either deformable 
or equipped with an elaborate cell adhesion model. 

In contrast to those, we are interested in a minimal 
resolution of RBCs since reducing the resolution gener- 
ally is the most effective way to enhance the efficiency 
of a fluid dynamics solver. Ahlrichs and Diinweg imple- 
ment a dissipative coupling of point-particles to an LB 
fluid |16|. However, the description of RBCs as point par- 
ticles would involve a resolution which is so low that hy- 
drodynamics in the smallest vessels becomes inaccessible. 
Concerning hydrodynamics, it is questionable whether 
the resolution of cell deformation has a benefit compared 
to a rigid particle model if each RBC is resolved with 
only a few lattice spacings. In consequence we decide for 
a method for suspensions of rigid particles of finite size 
that is similar to the one described in [13| . As will be ex- 
plained below, not volume-conserving cell deformations 
of the order of one lattice spacing occur as an artefact 
of the method already but do not show significant infiu- 
ence on the fiow behavior. Ding and Aidun simulated 
rigid particles with the biconcave shape of unstressed 
red blood cells using an LB method [1^. It is known, 
however, that RBCs abandon this equilibrium shape and 
instead resemble elongated ellipsoids when exposed to 



shear flow [l9|. Thus, taking into account the limited lat- 
tice resolution, we decide for discretized ellipsoidal model 
cells with rotational symmetry as an approximation of 
the shapes actually assumed by real eryth rocy tes in many 
flow situations. Differently from [H, |20, l21|, our imple- 
mentation does not enforce rigid particle surfaces since 
this would be in contradiction with the nature of de- 
formable erythrocytes. Especially in bulk flow at high 
volume concentrations but also in capillaries due to the 
influence of walls, the flow is not dominated by long-range 
hydrodynamics but by short-range cell-cell and cell-wall 
effects. Thus, a coarse-grained description using effec- 
tive cell and wall interactions is appropriate. We account 
for the complex short-range behavior of RBCs on a phe- 
nomenological level by means of model potentials. Our 
potentials serve to provide a softly repulsive core that fol- 
lows the approximated ellipsoidal RBC shape. For this 
purpose, the method of Berne and Pechukas [22| is ap- 
plied in order to anisotropically rescale a Hookian spring 
potential. 

In the following section |ll] we provide an introduction 
to the application of the LB method to suspensions of 
rigid particles and discuss how our model differs from 
other implementations. In section IIIII we develop phe- 
nomenological potentials for the anisotropic interaction 
of two cells and of cells and walls. Section HVl opens with 
the search for a parametrization which fits to experimen- 
tal literature data. This set of parameters is then used to 
demonstrate the applicability of the new model to con- 
fined systems. We further discuss the performance of 
our implementation for large systems and conclude with 
a summary in section [V] 



II. HYDRODYNAMIC PART OF THE MODEL 

To model the blood plasma that surrounds the cells the 
LB method is applied. For a comprehensive introduction 
we refer to [lOJ . The fiuid traveling with one of r discrete 
velocities c^ at the three-dimensional lattice position x 
and discrete time t is resembled by the single particle 
distribution function ^^(x, i). Its evolution in time is 
determined by collision 



ri*(x, t) = nr{x., t) — fl 
and the successive advection 

nr{y. + Cr,t + 1) = n,*(x, t) 



(1) 



(2) 



of the post-collision distribution n,*(x, t). Eq.[T]and Eq.[5] 
together can be written as the lattice Boltzmann equa- 
tion 



nr(x -|- Cr, t + 1) = ?lr(x, t) — fl . 



(3) 



For the sake of siinplicity and computational efficiency we 
follow a D3Q19 approach with a single relaxation time 



r [231 ■ We thus have 19 discrete velocities and the BGK- 
colhsion term [2J| 

^ ^ n^(x,0-<q(p(x,t),u(x,t)) 



To model boundaries moving with velocity v, Eq. [5] 
can be complemented with a correction term 
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where the equilibrium distribution function 



n';'^{p,u) = pac 



CrU (CrU) U^ 



(CrU) U^CrU 



6c6 



2c4 



(5) 



is an expansion of the Maxwell-Boltzmann distribution 
of third order in velocity u [2^ . The value of the speed 
of sound Cs = 1/a/3 depends on the choice of the lattice. 
The same holds for the lattice weights 



(6) 



which differ for lattice velocities c^ according to their 
absolute value Cr- The local density 
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P{^,t) =^rv(x,i) 



and velocity 



u(x,t) 



J2rnr{'^,t)Cr 

p(x,i) 



(7) 



(8) 



are calculated as moments of the fluid distribution with 
respect to the set of discrete velocities. Both are invari- 
ants of the BGK collision rule Eq. [TJ This method is well 
established for the simulation of the liquid phase of sus- 
pensions [13; [201 1 namely blood [y, Q. It can be shown 
that in the limit of small velocities and lattice spacings 
the Navier-Stokes equations are recovered with a kine- 
matic viscosity i' = (2r — l)/6. 

For a coarse-grained description of the hydrodynamic 
interaction of cells and blood plasma, a method similar to 
the ones explained in [17| and [20| modeling rigid parti- 
cles of finite size is applied. Starting point is the mid-link 
bounce-back boundary condition that implements no-slip 
boundaries for the fluid; arbitrarily shaped geometries 
are discretized on the lattice by turning the lattice nodes 
on the solid side of the theoretical solid-fluid interface 
into fluid-less wall nodes. If x is such a node the up- 
dated distribution at x -I- c^ is not determined by the 
advection rule Eq. [2] but according to 



nr{x + Cr,t +1) = n^(x + Cr,t) 



(9) 



which means replacing the local distribution with the 
post-collision distribution of the opposite direction f (de- 
fined by Cf = — Cj.). We make use of this boundary con- 
dition to implement (rigid) vessel walls. 



which is of first order in velocity. Inserting Eq. [5] and 
Eq. [ini one can easily prove that the new update rule 



nr{x + Cr,t+ 1) = n,*(x + Cr,t) + C 



(11) 



is up to second order consistent with the equilibrium dis- 
tribution function Eq. [S]for the general case u = v 7^ 0. 
When used to implement freely moving particles, it is 
necessary to keep track of the momentum 



^Pfp = (2"-f + G) Cr 



(12) 



which is transferred during each time step by each single 
bounce-back process. According to the choice of unit 
time steps it is equal to the resulting force on the particle. 
The equations of motion of the particles are integrated 
like in classical molecular dynamics (MD) simulations to 
achieve the time evolution of the system. We implement 
a combined LB/MD code in which both the time step and 
the spatial decomposition scheme are shared between the 
two methods. A leap-frog integrator that is adapted to 
the internal representation of the cell orientations based 
on unit quaternions is applied |26|. 

Due to discretization errors the representation of a par- 
ticle on the lattice slightly changes its shape and volume 
during movement. When new lattice sites are covered, 
the fluid at those is deleted. When a site formerly occu- 
pied by a particle is freed, new fluid is created according 
to Eq. [3 In doing so, the initial density p of the simula- 
tion is utilized as p. The velocity u is estimated according 
to the translational and rotational velocity of the particle 
and the no-slip assumption. In both cases the change in 
total fluid momentum is balanced by an additional force 
on the respective particle. 

Physiological RBCs, however, are deformable and as- 
sume the shape of biconcave discs in the absence of exter- 
nal stresses [^. Despite the coarse-graining in our model, 
we do not want to give up the anisotropy of RBCs. Obvi- 
ously, anisotropic model cells are able to display a much 
richer behavior than radially symmetric particles. We 
thus choose a simplified ellipsoidal geometry that is de- 
fined by two distinct half-axes R\\ and R± parallel and 
perpendicular to the unit vector 6^ which points along 
the direction of the axis of rotational symmetry of each 
particle i. 

Closely approaching particles are modeled as follows: 
when there are still fluid nodes between both discretized 
volumes the LB method is able to keep track of the 
emerging lubrication forces apart from discretization er- 
rors. As soon as there is a direct particle-particle in- 
terface without intermediate fluid nodes, the lubrication 
forces cannot be covered by a lattice-based method any- 
more. Moreover, an effective attraction becomes visible 



because of the missing fluid pressure in between the par- 
ticles. Typical applications of this simulation method 
to the case of dense suspensions additionally feature an- 
alytical short-ran ge l ubrication corrections to overcome 
this problem [20l . l2l| . These are implemented as pair- 
forces that depend on the relative velocity and diverge 
for vanishing gap-widths. However, this procedure is in- 
appropriate for a model for suspensions of deformable 
cells. Since the theoretical particle shapes defined by i?|| 
and R± are fixed, our application even requires toler- 
ance for the overlap of the discretized volumes in order 
to account for the case where two cells strongly deform 
while approaching each other. Due to the complexity 
of the emerging forces that include electrostatic repul- 
sion and van der Waals forces but also the mechanics 
and chemistry of the cell membranes and the rheology of 
the cell plasma, we cover them on a purely phenomeno- 
logical level in the next section. Here, we support the 
LB method with additional rules which result in forces 
for the case of two particles in direct contact with each 
other that are neither divergent nor excessively attrac- 
tive. Wherever a direct particle-particle interface is en- 
countered, we apply a pair of mutual forces 



and 



F+p = 2<<i(p,u==0)c, 



Fpp = 2nt\p, u = 0)cp = -F 



pp 



(13) 



(14) 



at each link across the interface which are directed to- 
wards each respective particle. Comparison with Eq. 1121 
shows that this is exactly the momentum transfer during 
one time step due to the rigid-particle algorithm for a 
resting particle and an adjacent site with resting fluid at 
equilibrium and initial density p. The fluid in our simu- 
lation is to good approximation incompressible and the 
velocities are small. The forces arising from those re- 
gions of the particle surfaces that are in contact with the 
fluid therefore are largely compensated and do not cause 
an artificial attraction. In consequence, the self-induced 
collapse of particles in contact is prevented without the 
need for divergent lubrication corrections as in rigid- 
particle models. Moreover, for a given system, Eq. 1131 
and Eq. [HI depend only on c^ = — Cp. For symmetry 
reasons, n^'^{p,0) = nf^{p,0) holds. Thus, the momen- 
tum balance is kept since the two forces emerging from 
any particle-particle link compensate each other. How- 
ever, since Eq.[T51and Eq. [T3]do not depend on the rela- 
tive velocity they cannot cover dissipative forces between 
particles. This limitation needs to be kept in mind when 
deciding about phenomenological cell-cell forces and their 
parametrization later in this paper. 

For the sake of simplicity we do not allow a lattice 
node to be occupied by more than one cell. Occupation 
instead is determined by the order in which cells arrive 
at a node. From the point of view of the surrounding 
fluid this behavior is physically consistent with two par- 
ticles that compressibly deform upon contact. The com- 
pressibility, however, can lead to an artificial increase of 



the total mass since the number of fluid nodes increases 
temporarily and we do not adjust the particle mass dy- 
namically according to the volume occupied momentarily. 
Thus, in an ensemble of many cells, there are fluctuations 
of the total mass due to the introduction of the correction 
term C in Eq. [TT] and due to the change in total volume 
of the solid phase which fluctuates when cells move and 
increases when they overlap. However, even during mil- 
lions of time steps we find no drift of the total mass for 
the systems we simulate here. 

In case of close contact of cells with the confining ge- 
ometry we proceed in a similar manner as for two cells. 
The only difference is that the forces on the system walls 
are ignored since they are assumed to be rigid and fixed. 



III. MODEL POTENTIALS FOR CELL-CELL 
AND CELL- WALL INTERACTIONS 



In order to account for the complex behavior of real 
RBCs at small distances we add phenomenological pair 
potentials. For simplicity, we restrict ourselves to repul- 
sive forces. This can be justified because in many phys- 
iological situations of interest, for example close to the 
walls of large parts of the arterial system, high shear rates 
render aggregative effects negligible [^, [23| . The task of 
the potential therefore lies in establishing an excluded 
volume for each cell. Due to the mild increase of the 
potential, an overlap of these volumes will be unfavor- 
able yet possible to some degree. Thus, the deformation 
of cells upon contact is modeled in a phenomenological 
way. A simplified potential also is beneficial to the effi- 
ciency of the model since it can be evaluated with less 
numerical effort and is less likely to demand small time 
steps or high order integrators. We start with the repul- 
sive branch of a Hookian spring potential 



e(l - nj/cr) 
"''~ ^ 






(15) 



for the scalar displacement r^j of two particles i and j. 
This is probably the simplest way to describe (elastic) 
deformability. The energy at zero displacement and the 
distance at which the repulsive potential force sets in can 
be directly controlled by means of the parameters e and 
a. With respect to the disc-like shape of RBCs, we follow 
the approach of Berne and Pechukas [22| and choose the 
stiffness parameter 



£(6„6j) 



and the size parameter 



1 - x^ (6,6j)' 



(16) 



a[6i,Oj,Vij) 



1 - ^ 

^ 2 






1-X6i6j 



(17) 



as functions of the orientations 6^ and 6j of the cells and 
their normalized center displacement f^. We achieve an 
anisotropic potential with a zero-energy surface that is 
approximately that of ellipsoidal discs. Their half-axes 
CTii and a±_ parallel and perpendicular to the symmetry 
axis enter Eq. \W\ and Eq. [17] via 



a ~ la ^ and x 



(18) 



II 



whereas e determines the potential strength. The above 
approach for anisotropic rescaling of radial symmetric po- 
tentials and its later improvement by Gay and Berne |28| 
were intended for modeling liquid crystal systems. Par- 
ticularly the method by Gay and Berne is applied almost 
exclusively to a Lennard- Jones potential featuring a short 
range repulsion and an attraction on moderate distances. 
This is referred to as "Gay-Berne potential" in the liter- 
ature. The model potential presented by us lacks the at- 
tractive tail but is equipped only with a softly repulsive 
core. In consequence, there is no force acting on parti- 
cles separated by more than the respective core diame- 
ter and at physiological volume concentrations we cannot 
expect to observe spontaneous ordering of the system. 
Compared to typical liquid crystal applications, the role 
of our potential lies rather in providing a soft repulsion 
within an anisotropic discoid volume than in making spe- 
cific cell alignments more favorable compared to others. 
Fig. [1] displays in dimensionless form the magnitude of 
the resultant repulsive pair force _F as a function of ry- for 
o\ 
6, 




FIG. 1: Dimensionless repulsive potential force as a function 
of the dimensionless center distance for a\ — Scrii and three 



sets of relative orientations 6i 11 6, 1. iij, 6i 



and 



6i -L 6j II Vij. An approximately ellipsoidal excluded volume 
can be deducted from the surface at which the repulsion sets 



instead of Eq. [T7] as a size parameter with 

; ^f - 

[ + al and Xw = ^t" 

"^11 " 



(20) 



3(T|| and three simple sets of relative orientations: 



^3 -'- ^ij I 



allows to scale a potential with radial symmetry to fit 
for the description of the interaction of a sphere and an 
ellipsoidal disc J22| . fix is the normalized center displace- 
ment of particle i and a wall node x. It is not necessary 
to scale the stiffness parameter anisotropically, instead 



and Oi -L Oj \\ Vij. Dependmg ^g gg^ e(6,;, 6^) = £w fixed and use Ew to tune the poten- 



on the orientations, the repulsive force sets in at differ- 
ent r.ij . Aiming at the presentation of a model potential 
which is simplified to the greatest possible extent, we 
chose the Berne-Pechukas approach which is slightly less 
complex than the more popular one by Gay and Berne. In 
this approach, it is not possible to independently adjust 
the interaction strength for different molecular orienta- 
tions. That the potential is considerably stiffer in the 
case where the flat sides of both cells are aligned towards 
each other is, however, consistent with the fact that for 
this orientation, the same linear approach creates a sig- 
nificantly larger overlap volume than in the other two 
cases. In the following section, we will find values for 
e, (Tji , and aj_ that reproduce the rheological behavior of 
blood. 

For modeling the cell-wall interaction we assume a 
sphere with radius a^ ~ 1/2 at every lattice node on 
the surface of a vessel wall and implement similar po- 
tential forces as for the cell-cell interaction based on the 
repulsive spring potential Eq. 1151 Berne and Pechukas 
show that using 



Cr 0,;,r,;x = 



1 - Xw (fixOi) 



(19) 



tial strength. The values of (T|| and a± are kept the same 
as for the cell-cell interaction. 

Fig. [2] shows a conclusive outline of the model. Two 
cells i and j surrounded by blood plasma and a section of 
a vessel wall are displayed. To enhance the explanatory 
power of the drawing we choose to restrict ourselves to 
the presentation of a cut parallel to the axes of rotational 
symmetry of the cells. Thus, the RBCs are visualized as 
two-dimensional ellipses instead of three-dimensional el- 
lipsoids. Depicted are the cell shapes defined by the zero- 
energy surface of the cell-cell potential Eq. [T5]with Eq. [T7] 
that can be approximated by ellipsoids with the size pa- 
rameters CT|| and cr_L as half axes. Also shown are the 
spheres with radius CTw defined accordingly by the cell- 
wall interaction Eq. [TS] and Eq. ^^ which are assumed at 
all wall nodes that are linked to a fluid node by one of the 
lattice directions c^. While these spheres are centered on 
the respective wall nodes, the cells are free to assume 
continuous positions and orientations o^ and Oj . In con- 
sequence, also the center displacement vectors r^ and r.^x 
between the cells and between cell i and an arbitrary wall 
node X are continuous. Still, for the cell-plasma interac- 
tion an ellipsoidal volume with half axes Ru and R± is 
discretized on the underlying lattice. For clarity, this is 
drawn only for cell j. 
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FIG. 2: (Color online) Outline of our 3D model by means 
of a 2D cut. Shown are two cells i and j with their axes 
of rotational symmetry Oi and Oj. The volumes defined by 
the cell-cell interaction is approximately ellipsoidal with half 



axes (T|j and a± (red, 



The ellipsoidal volume of the cell- 
plasma interaction with half axes _R|| and R± is discretized on 
the underlying lattice. It is shown for only one cell (blue, - -). 
The cell-wall potential assumes spheres with radius CTw on all 
surface wall nodes (green, — ). Depicted are also the center 
displacement vectors r^j and r^x between both cells and to an 
arbitrary surface wall node x. 



IV. RESULTS 

As a convention in this work, primed variables are 
used to distinguish quantities given in physical units from 
the corresponding unprimed variable measured in lat- 
tice units. The maximum extent of physiological RBCs 
at their equilibrium shape amounts to about 8 ^m and 
2.6 ^m perpendicular and parallel to the axis of rotational 
symmetry [2| ■ We find that an ellipsoid of revolution with 
the same numbers as axes has a volume of 87 ^m'^ which 
fits with the RBC volume measured in [2|. We therefore 
choose the size parameters of the cell-cell potential to be 



(7^ = 4/ini and ctm = 4/3 /im 



(21) 



and achieve that both the magnitude and the maximum 
extents of the volume defined by the cell-cell interaction 
match typical values for physiological erythrocytes. 

All quantities that are of interest in our simulations 
can be converted from simulation units to physical units 
by multiplication with products of iirteger powers of the 
coirversion factors Sx, St, and Sm for space, time, aird 
mass that thus completely define a scale. We determine 
the mean deviations of the Stokes drag coefficients of a 
single spherical particle from the theoretical values in the 
laminar regime to be in the order of 10~^ for a radius 
of 2.5 lattice units. This is in agreement with equivalent 
tests done by Ladd [29| . We therefore restrict ourselves to 
simulations of particles whose representation on the lat- 
tice is at least as large as that of a sphere with radius 2.5. 
When using the same aspect ratio R^/R\\ = ai_/a\\ = 3 



for the cell-fluid as for the cell-cell interaction, this re- 
quirement results in minimum values for R± and _R|| of 
3.6 and 1.2 lattice units, respectively. 

It can be expected that with cell-fiuid volumes that 
are significantly smaller than the size parameters of the 
potential we cannot achieve realistic coupling strengths 
which are needed for example to model the clogging of 
capillaries. Still, R± and Ru should be smaller than 
the respective size parameters of the cell-cell potential 
since limiting the amount of overlapping cell-fiuid inter- 
action volume will improve the modeling of hydrodynam- 
ics between cells. Ladd et al. |30| suggest assisting the 
particle-fluid coupling method with lubrication correc- 
tions starting at gap widths below 2/3 lattice spacings. 
Throughout this work, we choose Sx = 2/3 /im as a com- 
promise that both keeps the resolution and the computa- 
tional cost low and allows one to combine — for example — 
a high ratio of -R||/cr|| ~ 7/8 with a minimum gap width 
of 2(cr|| ~ R\\) — 0.5 at which the cell-cell potential starts 
to set in. 

To improve the numerical stability of the LB method 
and to easily relate given input radii i?|| and R± to an 
effective particle size [23, we always set the relaxation 
time to T = 1. This, together with the constraint 
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caused by the fact that the simulated kinematic fluid 
viscosity ly is supposed to match the kinematic viscos- 
ity of blood plasma of i^' = 1.09 x 10~^m^/s when con- 
verted to physical units, determines the time discretiza- 
tion as St = 6.80 X 10~* s. For convenience, we arbi- 
trarily choose the fluid density in simulation units to 
he p ~ 1. With Sx and the physical plasma density 
p ~ 1.03 X lO^kg/m'', this choice results in a mass con- 
version factor of Sm — 3.05 x 10~^^ kg. 

We first investigate the effects of the free model pa- 
rameters by measuring the ratio of the apparent dynamic 
viscosity /iapp and the constant plasma viscosity /i for a 
homogeneous suspension of cells in plane Couette flow. 
All simulations reported here are performed on a sys- 
tem with a size of l^, = 128 lattice units in x- and at 
least ly ~ Iz = 40 lattice units in y- and z-direction. 
This represents 85 x 27^ /im^ of real blood. Between the 
two j/2-side planes a constant offset of the local fluid 
velocities in z-direction is imposed by an adaption of 
the Lees-Edwards shear boundary condition to the LB 
method [3l|, [33 . The other edges are linked purely peri- 
odically. For the cells, we implement a reflective bound- 
ary condition that negates the normal velocity compo- 
nent of a cell as soon as its center distance to one of the 
sheared side planes becomes less than a±_ . This procedure 
surely is inconsistent with respect to the open boundaries 
implemented for the fluid but far easier to achieve than a 
common Lees-Edwards implementation for both phases. 
To prevent these boundaries from influencing our mea- 
surements, we determine the shear rate 7 only in the 
central half of the system where the flow resembles an 



unbounded Couette flow. We obtain 7 from a linear fit 
of the velocity profile v^ {x) . The apparent viscosity 
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is then calculated based on A^le.z which is the averaged 
z-momentum transfer across both Lees-Edwards bound- 
aries during one time step. For each shear rate, we 
start with resting and randomly oriented model cells sus- 
pended in likewise resting fluid. We calculate Eq. [23] in 
intervals during the simulation and start accumulating 
the result for temporal averaging as soon as a steady 
state is achieved. Several samples prove that neither the 
change of the random seed for the generation of the initial 
cell configuration nor the stepwise increase of the system 
size perpendicular to the velocity gradient up to a vol- 
ume of 85'^ /im"^ leads to any significant deviation of the 
results. However, we find that the shear causes the cell 
orientations {6^} to preferably align in the xz-plane. 

A proper choice of the ratio Rn/an is not known a pri- 
ori. We thus perform simulations at a constant shear 
rate of 7 = (2.21 ± 0.08) x 10^ s^^ for different i?j|/cr||. 
The resulting particle Reynolds numbers i?ep are of the 
order of 10~^. We arbitrarily choose a cell number den- 
sity of p' = (6.4 ± 0.3) X 10^^ m^'^ corresponding to a 
physiological hematocrit of 56 % and a cell stiffness pa- 
rameter of e' = 1.47 X 10~^^ J. The resulting apparent 
viscosity //app as a function of the ratio R\\/a\\ is drawn 
in Fig. [3l A relatively mild and almost linear increase 
is visible for R\\ /a\\ < 1 which can be related to the in- 
crease of friction in the system. Around 1 . the increase 
becomes considerably steeper as the minimum gaps of 
approaching cells vanish. At even larger ratios, the slope 
decreases again due to large and unphysical amounts of 
overlap of the cell-fluid coupling volumes that accordingly 
to Eq . [Ol and Eq. [Til reduce the effective friction between 
cells. Based on our previous considerations and affirmed 
by Fig. [31 we choose R\\/<Ji\ = 11/12 sa 0.92 as a value 
that is close to unity but still induces an only moderate 
amount of overlap even at high shear rates of the order 
of 10'^ s~^. This choice results in size parameters of the 
cell-fluid interaction of 



i?l = 11/3 ,um and R\, = 11/9 ^tm 



(24) 



All dimensions in Fig. [5| above were already drawn to 
scale with respect to the dimensional parameters in 
Eq. [H] and Eq. [Mj 

The parameter s is of special interest since it con- 
trols the cell stiffness which describes the deformabil- 
ity of the erythrocytes in our model. From experiments 
it is known that the shear thinning behavior of blood 
at high shear rates is related to the deformability of 
the RBC membrane and can be disabled by artificial 
hardening of the cells [3, [Sj]- Our implementation of 
the model stays numerically stable only for a limited 
range of e. Simulations performed for various shear rates 
1.7 X 10^ s~^ < 7 < 2.3 X 10^ s~^ corresponding to par- 
ticle Reynolds numbers 10~^ < Rep < 1 and e' varying 
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FIG. 3: Dependence of the apparent dynamic viscosity /iapp 
on the fraction R\\/<y\\ of the linear dimensions of the cell- 
fluid and cell-cell interaction volumes for a shear rate of 7' = 
(2.21 ±0.08) X 10^s"^ a number density of p' = (6.4 ±0.3) x 
10^^ m~^, a cell stifl^ness parameter e' = 1.47 x 10~^^ J, and 
cell-cell size parameters a'^ — 4 ^m and a'« — 4/3 /xm. All 
consecutive simulations are performed with -R||/o"|| ~ 11/12 « 
0.92. 



between 1.47 x IQ-i^^ J and 1.47 x 10"^^ j ^^ ^ cell-fluid 
volume concentration of 43 % still show that for a given 
shear rate, larger e result in higher viscosities yet in a 
less steep viscosity decrease. Fig. [3| displays this effect 
which is asymptotically consistent with the experimental 
results of Chien Q . It is interesting to note that by plot- 
ting the apparent viscosity over the fraction j/e — as we 
do in Fig. |4(b)| — we can collapse the region of strongest 
viscosity decrease in the curves for different e. This in- 
dicates that the shear thinning is determined by a bal- 
ance of viscous and potential forces that scale with 7 
and e, respectively. Comparison of Fig. [4| with experi- 
mental data taken from the literature [9| shows best con- 
sistency in the case of high shear rates 7' ~ 10"^ s~^ for 
e' = 1.47 X 10~^^ J. We keep this value for all further 
simulations in the current work. 

Having defined values for all parameters of the cell- 
fiuid and cell-wall interaction, we can now investigate 
the effect of varying cell concentrations on the viscosity. 
For given i?|[ and R±, the cell-fluid volume concentration 
<f>cf is proportional to the number concentration. Fig. [5| 
shows the dependence of the apparent viscosity on $cf 
for a fixed shear rate of 7' = (2.2 ± 0.1) x 10^ s'^. The 
particle Reynolds number is of the order of Rep ^ 10^^. 
For <I>cf < 35 % we find a nearly linear increase of //app- 
For <f>cf > 35%, the curve is still linear but the slope is 
slightly smaller. Compared to the literature, /iapp stays 
clearly below the viscosities known for hard ellipsoids 
with a sinfilar aspect ratio of 0.3 [3J|. The lower vis- 
cosities of our model, especially at high volume fractions, 
are caused by the reduced dissipation between touching 
and overlapping cells. This explanation can be substan- 
tiated by considering two neighboring boxes of a peri- 
odic arrangement each containing one cell in the center. 
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FIG. 4; (a) shows the shear rate dependent apparent viscosity 



/iapp at a cell- fluid volume concentration of $cf = 43 %. The 
different symbols represent different cell stiffness parameters 
e' = 1.47 X 10"'= with k = 16, 15, 14, 13, 12 from bottom to 
top. Rescaling the shear rate 7' with e' as displayed in (b) 
leads to a collapse of the region of steepest viscosity decrease 
on a single curve which hints at a concurrence of viscous and 
potential forces. All further simulations are performed with 
k = 15. 



Depending on the orientation and offset relative to the 
LB grid, direct cell-cell links start to occur at volume 
concentrations between 30 % and 50 % for the given R± 
and _R|| . These numbers also match the region in Fig. [S] 
where the slope of /Xapp(<i>cf) decreases. As can be seen 
from Fig. [SI our results for concentrations up to about 
50 % fit reasonably well with the experimental studies of 
Goldsmith (see [1|) and of Shin et al. [Hi. At higher $cf, 
touching cell-fluid volumes start to dominate the rheol- 
ogy of the model suspension. The exact shear rates ap- 
plied by Goldsmith and Shin et al. are not known. We 
can only infer from the literature that 7' was larger than 
100 s~^ and 250 s~^, respectively. In this range, blood 
shows shear thinning behavior [3, [S^l and so does our 
model (see Fig. |4]). It therefore is not possible to deter- 
mine whether — as Fig. [S] suggests — the model perfectly 
matches with experiments for $cf < 40 % and underes- 
timates the viscosity between 40% and 50%. However, 
Fig. 0] demonstrates that a better consistency at physi- 
ologically important concentrations around 40 % should 
be easily attainable by tuning the value of e. 

While the previous simulations regard bulk properties 
we now turn to examples where confinement and par- 
ticulate effects play a crucial role. The cell-wall inter- 
action stiffness e^ can be determined similarly as e by 
comparison with experimental data. As an example we 
choose the sieving experiments performed by Chien et 
al. who filtered human erythrocytes through polycarbon- 
ate sieves with mean pore diameters of D' — 2.2 /im 
to 4.4 /im and mean pore lengths of 13 /xm [33. They 
analyzed the resulting fiow resistance and damaging of 
cells in dependence on the pressure drop AP' across the 
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FIG. 5: fJ.app dependence on the volume concentration "l>cf 
related to the cell-fluid interaction at a fixed shear rate of 
(2.2 ± 0.1) X 10"^ s~^ compared to experimental data for 
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7' > 100 s^^ given in ^,^. 



sieves which was varied between approximately 10^ and 
lO^N/m^. We simulate a single cell in front of a pore 
at a smaU value of AP' = 4 x lO^N/m^ (0.3cmHg) and 
vary the pore diameter and s^. At this pressure drop, no 
significant hemolysis, which — as a sub-cell effect — is not 
resolved in our model, was found in the experiments |35| . 
However, compared to the case of D' — 4.4 ^m, the flow 
resistance was increased by a factor of approximately 4 
for D' = 3.7 /im, by about 30 for D' = 3.0 /im, and by 
more than 100 for D' = 2.2 /im at a hematocrit not higher 
than of the order of 1 %. In our simulations, we iden- 
tify the non-passing of model cells with a high increase 
of flow resistance in the experiments. We find that for 
e'^ = 1.47 X 10~^^J, the cell passes a pore with only 
D' = 3.0 /im while for e'^ = 1.47 x 10"" J already a di- 
ameter of D' = 4.4 /im proves an insurmountable obsta- 
cle. In view of reference [351 these two Ew are unrealistic 



but the intermediate value e' 



1.47xl0~i'5jisan 



appropriate choice for this setup which allows the model 
cell to pass through pores with a diameter of 3.7 /im and 
more. We now study the fiow through a bifurcation of 
cylindrical capillaries with a radius of 4.7 /xm. One of 
the branches, however, features a stenosis with radius 
Rs- A cut through the geometry containing nine RBCs is 
displayed in Fig. [Sj It visualizes the cells as the approxi- 
mated ellipsoidal volumes defined by the zero-energy sur- 
face of the cell-cell interaction and the vessel walls as 
midplane between fiuid and wall nodes. The open ends 
of the system are linked periodically. We drive the sys- 
tem by means of a body force acting only on the fluid in 
the entrance region. As initial condition, cells are placed 
randomly in the unconstricted parts of the system. Both, 
the tube diameters and Reynolds numbers Re ^ 4 x 10"'^ 
match physiological situations Q. As above, the cell- wall 
potential stiffness is chosen to be e^ = e"' = 1.47x10"^^ J. 
We average the relative flow rate through the constricted 

branch yconstr = yconsti/ (Qconstr + tt/unconstr) irOm 1.7 S 
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FIG. 6: Cut through a capillary bifurcation. Shown are the 
volumes defined by the cell-cell interaction and the midplane 
between wall and fluid nodes. The plasma is not visualized. 
The flow direction is fi'om left to right. The vessel radius is 
i?s = 2.7 /xm at a stenosis in the upper branch and 4.7 /^m 
otherwise. Geometries with length scales that are not large 
compared to a cell diameter require treatment by a method 
that is able to resolve particulate effects like the shown clog- 
ging of the constricted branch for a cell- wall potential strength 
of e' = 1.47 X 10"^^ J. 



0.5 


- -H 


total ' ' 


. i 


_ 






plasma A: ' 


d) 


- 


■i 0.4 

CD 


cells y\--Y 


J. 


1 0.3 
i 0.2 


~ 


j/^ 




~ 


relati 
o 


■y 


^ 




_ 




9 


Q 


□ 


- 


0.0 


- B 


1 1 1 


■" 



2.5 



3.0 



3.5 4.0 



4.5 



5.0 



FIG. 7; Time-averaged relative flow rate through the con- 
stricted capillary in a bifurcation as shown in Fig. [6] for dif- 
ferent stenosis radii _Rs. The cell- wall interaction stiffness is 
e^ = 1.47 X 10"^^ J. While for i?^ = 2.7 ^m, the constricted 
branch becomes clogged and only a small amount of plasma 
passes the remaining aperture, no clogging occurs for larger 



to 3.0 s measured from system initialization. This is done 
for i?g = 2.7 /xm, 4.0 /im, and 4.7 /im. As expected, the 
results in Fig.[7]are monotonous with R^. When studying 
the volumetric flow rates of plasma and cells separately, it 
becomes clear that for E!^ = 2.7 /Ltm the cells cannot pass 
the constriction at the present body force. This situation 
is visualized in Fig. [51 

We further study the dynamics of the system for the 
present and two lower cell-wall interaction parameters 
and plot QconstrW in Fig. ^ For ^ = 1.47 x IQ-i^ J, 
the curve decreases in two steps due to the successive 
arrival of two erythrocytes and stays below 10 % as only 




FIG. 8: Time evolution of the relative flow rate through the 
constricted capillary in the bifurcation shown in Fig. [6l For 
a cell-wall interaction stiffness of E'„ = 1.47 x 10^^^ J, the 
relative flow rate drops to less than 10 % in two major steps 
related to the successive arrival of two single erythrocytes at 
the constriction. With e^ — 1.47 x 10^^® J, the reduction of 
the flow rate is only temporary, since the cells are eventually 
able to pass. While leaving the constriction, the RBCs are 
accelerated by the cell-wall potential forces which explains 
the peaks of the flow rate. 



a small amount of plasma is able to pass the remain- 
ing aperture. In contrast, there is a continuous flow of 
plasma and cells for e'^ — 1.47 x 10~^^ J. For an inter- 
mediate stiffness e^ = 1.47 x 10"-'^^ J the cells get stuck 
initially. However, the flow in the narrowed branch is in- 
fluenced by the time-dependent cell configuration in the 
other branch. It happens eventually that the pressure in 
front of the stenosis rises to a level which lets the RBC 
overcome the barrier imposed by the cell-wall potential. 
Each restitution of a higher flow level is initiated by a 
peak which can be explained by the cell-wall potential 
that accelerates the RBC into the flow direction while 
the cell leaves the constriction. As another effect, we 
flnd Ew to affect the relative flow rates of the two phases 
since larger values force the cells into the center of the 
capillaries where higher velocities are measured. 

Despite the coarse-graining of the model it qualita- 
tively reproduces some aspects of the behavior observed 
for blood flow in capillaries. The fact that cells approach- 
ing a bifurcation show a strong preference to choose the 
faster branch is described in [3[. The last consequence 
of this effect is visible in Fig. [6] and Fig. [8] where dur- 
ing a considerable amount of time no further RBCs en- 
ter the constricted branch after its closure. It can be 
seen that our model is able to describe particulate effects 
which could hardly be covered in terms of a continuous 
fluid. Obviously, reproducing the behavior of single cells 
at bifurcations is crucial if the microcirculation and its 
heterogeneous flow properties are to be modeled [a|. 

The vessel radii present in the human microvascular 
system approximately cover a range from 2 /im to 50 [iva.. 
Having demonstrated the applicability of our model to 
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small capillaries, we proceed with a study of the steady 
flow through a larger vessel with a radius of i?' = 31 /im 
corresponding to an arteriole or venule [5| . In the simula- 
tion, the vessel is closed periodically at a length of 43 ^m. 
We choose <I>cf = 42 % and an intermediate cell- wall in- 
teraction stiffness e^ = 1.47 x 10^^^ J. Fig.[9]shows a cut 
through the vessel for steady flow at Re ^ 1. The flow 
is driven by a body force which acts on both plasma and 
cells in the whole system and is equivalent to a constant 
macroscopic pressure gradient. The system is evolved in 
time until neither the initial fee ordering of the cells nor 
significant directed changes in the volumetric flow rate 
Q are visible. In Fig. [TOl the radial velocity proflle in 
the case of a body force resulting in a pseudo-shear rate 
of v'^ = Q/(2itR^) = 1.3 X lO^s-i or a Reynolds num- 
ber of i?e '--^ 10 is shown. The graph deviates from the 
parabolic Hagen-Poiseuille profile that could be observed 
for a Newtonian fluid. Instead a parabolic core region 
and a narrow boundary region with high shear rates can 
be identified. The fit in Fig. [TU] shows that this pro- 
file can be easily explained by the modified axial-train 
model as described by Secomb [3a|- The fit parameters 
are the viscosity ratio of core and boundary /ic/Mb and 
the width of the cell-depleted boundary layer 5. The 
obtained viscosity ratio of fic/fJ-h ~ 2.43 ± 0.01 is con- 
sistent with the bulk properties in Fig. [5] if we assume 
fih = fJ. and 0.4 < $cf < 0.5 in the core. Also our result 
of S' = (1.47 ± 0.04) ^m seems compatible with the value 
of 1.8 ^m suggested by Secomb |36|. In additional stud- 
ies of radial cell-fluid concentration proflles we prove the 
existence of a cell-depleted layer and the possibility to 
tune its width by the cell- wall potential stiffness e^. We 
also flnd an increased cell concentration of up to around 
$cf = 60 % close to the central axis of the vessel. This 
must be a collective effect since in consistency with a 
2D study by Qi et al. [331 , we observe that single cells in 
Poiseuille flow migrate to an intermediate lateral position 
between vessel wall and center. 

The apparent viscosity for a cylindrical vessel is calcu- 
lated as 



TTR'^dP 



Mapp 



(25) 



with dP/dz being the macroscopic pressure gradient [5|. 
Pries et al. [3^ present an empirical expression for the 
dependency of /Xapp on the radius and hematocrit for the 
case of high flow rates after combining a variety of ex- 
perimental studies for pseudo-shear rates v'^ > 50 s~^ 
in a single flt. We perform a series of simulations at 
R' — 31 /im and three flxed pseudo-shear rates between 
v'^ = Q/{2ttR^) = (62±l)s-i and (563 ± 3) s"! but 
varying cell-fluid volume concentrations $cf- The corre- 
sponding Reynolds numbers are Re ^1. If we identify 
$ct with the hematocrit as we do in Fig. [5J we find very 
good agreement with the relationship by Pries et al. (Sg ■ 
The comparison is plotted in Fig. 1111 

The presence of a cell-depleted layer is closely con- 
nected to the emergence of heterogeneous cell concen- 
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FIG. 9: Cut through a cylindrical vessel with a radius of 
31 /im. For this geometry we choose a cell-wall interaction 
strength of e(v — 1.47 x 10~^^ J. Shown are the volumes de- 
fined by the cell-cell interaction at 42 % cell-fluid volume con- 
centration and the midplane between wall and fluid nodes. 
The flow is pointing into the drawing plane and has a maxi- 
mum velocity of 1.08 x 10^^ m/s at the center. 
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FIG. 10: Radial velocity profile in a cylindrical vessel with 
42 % average cell-fiuid volume concentration. Apparent slip 
due to a cell depletion layer is visible. The profile can be 
well fitted by a modified axial-train model as described by 
Secomb [3q |. The parabolic Hagen-Poiseuille profile is plotted 
as well for comparison. 



trations in different parts of the microvasculature since 
branching daughter-vessels first of all drain blood from 
the boundary layer [3|. The hematocrit, in turn, influ- 
ences the flow resistance, the flow rate, and the resulting 
distribution of erythrocytes at branching points [5J . Our 
simulation approach reproduces these aspects at least 
qualitatively. When implemented together with an in- 
dexed LB scheme as in |11[fL3| . the method will be able 
to simulate flow through digitized vessel networks cov- 
ering the whole scale of the microcirculation with high 
efficiency. Such simulations are still computationally de- 
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FIG. 11: Dependence of the apparent dynamic viscosity /iapp 
in a cylindrical vessel with radius B! = 31 /im on the cell-fluid 
volume concentration $cf- Three pseudo-shear rates v'^ — 
Q/(2ttR^) are examined. The empirical result by Pries et 
al. |38l | for v'^ > 50 s~^ with "l?cf as tube hematocrit is plotted 
for comparison. 



manding despite the simplifications of the model. Thus 
even systems that are small in physical units require par- 
allel supercomputers which makes the scalability of the 
code crucial. For a quasi-homogeneous chunk of suspen- 
sion consisting of 1024^ x 2048 lattice sites and 4.1 x 10^ 
cells simulated on a BlueGene/P system, we achieve a 
parallel efficiency normalized to the case of 2048-fold par- 
allelism of 95.7% on 16384 and still 85.2% on 32768 
cores. In comparison, the pure LB code without the 
MD routines that are responsible for the description of 
the cells shows a relative parallel efficiency of 98.1 % on 
32768 cores. The parallel performance of the combined 
code is mainly limited by the relation of the interaction 
range of a cell to the size of the computational domain 
dedicated to each task. We are aware of only one work 
on simulations of comparably large systems with a par- 
ticulate description of hemodynamics. This work was 
published by the group of Aidun and models the de- 
formation of cells explicitly [^, [3^. However, owing to 
the coarse-graining, our model is easier to parallelize ef- 
ficiently and — compared to the resolution given in [39| — 
allows for substantially higher cell numbers. Generally, 
our relatively low spatial resolution is highly beneficial for 
the simulation of large systems since from Eq. [51] it can 
be derived that the number of lattice site updates neces- 
sary for the simulation of a system with a given physical 
size for a given physical time interval scales with the fifth 
power of l/5x. As for plain LB simulations, this number 
is a good measure for the computational cost also in the 
case of our suspension model. 



V. CONCLUSION 

We developed a new approach for the coarse-grained 
simulation of suspensions of soft particles. This approach 
is based on a well-established method for rigid particle 
suspensions |17l |30| which covers the hydrodynamic long- 
range interactions and phenomenological model poten- 
tials to account for the behavior at small particle sep- 
arations. A parametrization suitable for the quantita- 
tive reproduction of hemorheology at moderate to high 
shear rates (j,[^,|33 was presented. The cell- wall interac- 
tion could be linked to experimental data on a single-cell 
level [33 ■ Afterwards, we demonstrated that the model 
shows a complex particulate behavior in bifurcations of 
partly constricted capillaries which is an essential fea- 
ture also of the flow properties of the microcirculation in 
vivo [5| . Using the example of steady flow through larger 
vessels, we proved the existence of a cell-depleted layer 
and obtained radial velocity profiles that are consistent 
with an accordant theoretical model [Sg • We could even 
quantitatively reproduce the experimentally observed de- 
pendency of th e ap parent viscosity in this geometry on 
the hematocrit [Sg- These results suggest that following 
our approach one can reproduce the particulate behav- 
ior of blood on a range of spatial scales that up to this 
moment was not covered by a single existing simulation 
method with comparable efficiency. Clearly, our motiva- 
tion is not to replace models with higher resolution like 
the ones presented in [a-|l; 113j but to bridge the gap 
to continuous descriptions of blood. We believe that this 
method can prove both an efficient tool for coarse grained 
yet particulate simulations of flow in microvascular vessel 
networks and a valuable contribution to the improvement 
of macroscopic blood modeling. 
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